1

Re-use your analysis and run an interaction model between your IV and the other covariate. Plot the interaction as a coefficient plot.
For the coefficient plot, you can stick to easystats.
linear_model_interaction <-
  lm(
    curfew_yes_no ~ age_cat * education_cat,
    data = gp_covid
  )

model_parameters(linear_model_interaction) %>% 
  plot()

2

Plot the interaction effect as a prediction.
For this purpose, you need the plot_model() function from sjPlot.
library(sjPlot)

plot_model(
  linear_model_interaction,
  type = "int"
)